function [rois,diffdata] = difference_segmentation(smalldata,params)


%% Automated segmentation of two-photon microscopy data
% moritz.kirschmann@fmi.ch 2013

%close all;
verbose=0; % 1 for results of median filtering, 2 for results of first segmentation


%% Loading of aligned dataset as smalldata

datasize=size(smalldata);
n_slices = datasize(3);


%% Definition of variables
if nargin ~= 2
    disp('using default parameters')
    meanoverdeltaT=5; % median filtering parameter, describing number of time steps
    
    sigma=1; % sigma of gaussan blur in xy
    
    deltaT=2; % delta t of the time derivative frames
    medianfilt = 0;
    medianfiltSize = 3;
    %set to 90
    vol_thres = 100; % volume threshold in pixel. All segments with less volume will be discarded
    %set to 90
    overlap = 0.5; % parameter of overlapping segements in 2d. If this overlap is reached the smaller object is discarded
    t_factor = 5;
    corrthres = 0.5; % threshold parameter of merging two segments with similar intensities
else
    disp('using custom parameters')
    meanoverdeltaT = params.meanoverdeltaT;
    sigma = params.sigma;
    deltaT = params.deltaT;
    vol_thres = params.vol_thres;
    overlap = params.overlap;
    corrthres = params.corrthres;
    medianfilt = params.medianfilt;
    medianfiltSize = params.medianfitlSize;
    t_factor = params.t_factor;
end 



%% Smoothing of raw data
disp('----smoothing----');
% mean filtering in time
ze = ones(1,1,meanoverdeltaT)/meanoverdeltaT; 

fdata=imfilter(smalldata,ze); 

% gaussian filtering in xy
xye = fspecial('gaussian',10, sigma);

ffdata=imfilter(fdata,xye);
%ffdata = fdata;
clear fdata;

%consider using imfilter on 3d with a 3d filter directly here
% h=fspecial3('gaussian',[5,5,5]);
% ffdata=imfilter(smalldata,h);
% 
meanIm=mean(ffdata,3);
meanIm=mat2gray(meanIm);

%% Creation of time derivativ data
disp('----calculating derivative----');
tic
diffdata =zeros(datasize(1), datasize(2), datasize(3),'double');

for i=meanoverdeltaT:n_slices-meanoverdeltaT
    diffdata(:,:,i) = ffdata(:,:,i+deltaT)-ffdata(:,:,i);
    if medianfilt == 1
        diffdata(:,:,i) = medfilt2(diffdata(:,:,i),[medianfiltSize medianfiltSize]);
    end
    
    if verbose == 1 
        figure(2);
        imagesc(diffdata(:,:,i), [-400 800]);
        figure(3);
        %imagesc(ffdata(:,:,i), clims);
        imagesc(ffdata(:,:,i))
        pause(0.5)
    end
 
end
toc

%% Oversegmentation of data
disp('----Oversegmentation----');
tic
absdiffdata=abs(diffdata);

me=mean(absdiffdata(:));
segmentation_threshold = t_factor*me;


mask = absdiffdata >= segmentation_threshold;
maskmat = reshape(mask,datasize(1),datasize(2),datasize(3));
toc


%% Removal of object with less pixels then vol_thres
disp('----Removing small segments----');
tic
nosmallmaskmat = bwareaopen(maskmat,vol_thres,6); 
disp('----Label unique regions / segments ----');
%[CC n_segs] = bwlabeln(nosmallmaskmat,6);
CC=bwconncomp(nosmallmaskmat,6);
n_segs = CC.NumObjects;

 if verbose == 2 
    for i=10: n_slices-10
        figure(20)
        imagesc(absdiffdata(:,:,i));
        figure(19)
        imagesc(maskmat(:,:,i));
        figure(22)
        imagesc(smalldata(:,:,i));
        
        figure(24)
        imagesc(CC(:,:,i))

        i
    end
 end
toc
clear absdiffdata;




%% Saving 2d masks of the segments in a cube for fast access, this

disp('----creating 2d maps from 3d data----');
tic
areaofmask = zeros(n_segs,1, 'uint16');
binxymask = false(datasize(1), datasize(2),n_segs);
 for n = 1: n_segs
     for j=1:numel(CC.PixelIdxList{n})
        [x,y,z]=ind2sub(CC.ImageSize,CC.PixelIdxList{n}(j));
        binxymask(x,y,n)=true;
     end
    
   

    areaofmask(n) = nnz(binxymask(:,:,n));
 end
 figure(22)
 mask=sum(binxymask,3);
 mask=bwperim(mask>0);
 overl=imoverlay(meanIm,mask,[0,0,1]);
 imshow(overl)
 toc
 pause(0.1)


 %% reduction of segments by looking which segments are contained
 %% within other segments.


 
 [ binxymask, n_redsegs ] = merging_overlap(n_segs,binxymask,overlap,areaofmask);
 
 figure(23)
 mask2=sum(binxymask,3);
 mask2=bwperim(mask2>0);
 overl=imoverlay(meanIm,mask2,[0,155,155]);
 imshow(overl)

 
 figure(24)
 overl2=imoverlay(overl,mask,[155,0,0]);
 imshow(overl2)
 title('r: before merging overlapping segments, b: after merging');
 toc
 pause(0.1)
 

rois = binxymask;
 end

